Models

In this document, we will fit all models regarding this project. The detailed preregistered analysis plan can be found in our OSF’s project.

Data

Let’s load libraries and data.

# Ensures the package "pacman" is installed
if (!require("pacman")) install.packages("pacman")

# Install/Load packages

pacman::p_load("rio",
               "here",
               "metafor",
               "dplyr",
               "bayesmeta",
               "brms")

remotes::install_github("stan-dev/cmdstanr")

# Ensure that packages are installed according to versions in renv's lockfile
renv::restore()

# Import data
d = rio::import(here::here("data", "data.xlsx"))

Let’s calculate the crude odds ratio based on the extracted data.

# Calculate crude log odds ratios

d_logOR = 
  metafor::escalc(
  measure = "OR", # log odds ratio,
  
  # Tocilizumab
  ai = trt_events,
  n1i = trt_total,
  
  # Control
  ci = control_events,
  n2i = control_total,
  
  data = d
) %>%
  as_tibble() %>% 
  # Calculate standard error
  dplyr::mutate(sei = sqrt(vi),
         # Set order to use "tocilizumab" as the Intercept in brms
         treatment = factor(treatment,
                         levels = c("tocilizumab",
                                    "sarilumab",
                                    "baricitinib")))

# saveRDS(d_logOR,
#         here("output", "data", "effect_sizes.Rds"))

Meta-analyses

First, we will fit a single meta-analysis model per drug (tocilizumab, sarilumab, or baricitinib vs. control).

All Bayesian meta-analysis random-effect models are the same, and is described as:

\[ \begin{align*} y_i & \sim Normal(\theta_i, \sigma_i^2) \tag{Likelihood} \\ \theta_i & \sim Normal(\mu, \tau^2)\\ \\ \mu & \sim \operatorname{Normal}(0, 0.75^2) \tag{Priors} \\ \tau & \sim \operatorname{Log-Normal}(-1.975, 0.67^2) \\ \end{align*} \]

where \(y_i\) is the observed mean log odds ratio in study \(i\) for either tocilizumab, sarilumab, or baricitinib versus control treatment. We assume these effect sizes are normally distributed around the study-specific mean \(\theta_i\) along with a known sampling variance, represented by the observed \(\sigma_i^2\). We also assume \(\theta_i\) is drawn from a normal distribution where \(\mu\) is the average effect and \(\tau^2\) is the between-study heterogeneity.

## Formula
mf = 
  # https://bookdown.org/content/4857/horoscopes-insights.html#consider-using-the-0-intercept-syntax
  formula(yi | se(sei) ~ 0 + Intercept + (1 | study))

## Priors

# Tau
informative = bayesmeta::TurnerEtAlPrior("all-cause mortality",
                                         "pharmacological",
                                         "placebo / control")


logmean = informative$parameters["tau", "mu"]
logsd = informative$parameters["tau", "sigma"]

# logmean
# -1.975

# logsd
# 0.67

priors_ma = 
  brms::prior(normal(0, 0.75), class = "b", coef = "Intercept") +
  brms::prior(lognormal(-1.975, 0.67), class = "sd") 

## Models

# Tocilizumab

ma_toci = 
  brms::brm(
  data = d_logOR |> dplyr::filter(treatment == "tocilizumab"),
  family = gaussian,
  
  formula = mf,
  prior = priors_ma,
  sample_prior = TRUE,
  
  control = list(adapt_delta = .95),
  backend = "cmdstanr", # faster
  cores = parallel::detectCores(),
  chains = 4,
  warmup = 2000, 
  iter = 4000, 
  seed = 123,
  
  file = here::here("output", "fits", "ma_toci.Rds"),
  file_refit = "on_change"
)

# Sarilumab

ma_sari = 
  brms::brm(
  data = d_logOR |> dplyr::filter(treatment == "sarilumab"),
  family = gaussian,
  
  formula = mf,
  prior = priors_ma,
  sample_prior = TRUE,
  
  control = list(adapt_delta = .95),
  backend = "cmdstanr", # faster
  cores = parallel::detectCores(),
  chains = 4,
  warmup = 2000, 
  iter = 4000, 
  seed = 123,
  
  file = here::here("output", "fits", "ma_sari.Rds"),
  file_refit = "on_change"
)

# Baricitinib

ma_bari = 
  brms::brm(
  data = d_logOR |> dplyr::filter(treatment == "baricitinib", study != "RECOVERY Bari"),
  family = gaussian,
  
  formula = mf,
  prior = priors_ma,
  sample_prior = TRUE,
  
  control = list(adapt_delta = .96),
  backend = "cmdstanr", # faster
  cores = parallel::detectCores(),
  chains = 4,
  warmup = 2000, 
  iter = 4000, 
  seed = 123,
  
  file = here::here("output", "fits", "ma_bari.Rds"),
  file_refit = "on_change"
)

Sensitivity analysis with all patients in RECOVERY Bari

ma_bari_sens = 
  brms::brm(
  data = d_logOR |> dplyr::filter(treatment == "baricitinib", study != "RECOVERY Bari (No Toci)"),
  family = gaussian,
  
  formula = mf,
  prior = priors_ma,
  sample_prior = TRUE,
  
  control = list(adapt_delta = .96),
  backend = "cmdstanr", # faster
  cores = parallel::detectCores(),
  chains = 4,
  warmup = 2000, 
  iter = 4000, 
  seed = 123,
  
  file = here::here("output", "fits", "ma_bari_sens.Rds"),
  file_refit = "on_change"
)

Meta-regressions

Now, we will fit Bayesian meta-regression models (one for tocilizumab vs. sarilumab data, another for tocilizumab vs. baricitinib).

The model can be described as:

\[ \begin{align*} y_i & \sim Normal(\theta_i, \sigma_i^2) \tag{Likelihood}\\ \theta_i & \sim Normal(\mu, \tau^2)\\ \mu &= \beta_0 + \beta_1 x_i\\ \\ \end{align*} \]

where \(y_i\) is the observed mean log odds ratio in study \(i\) for either tocilizumab, sarilumab, or baricitinib versus control treatment. We assume these effect sizes are normally distributed around the study-specific mean \(\theta_i\) along with a known sampling variance, represented by the observed \(\sigma_i^2\). We also assume \(\theta_i\) is drawn from a normal distribution where \(\mu\) is the average effect and \(\tau^2\) is the between-study heterogeneity. We will estimate \(\mu\) as a function of the estimated population effect size \(\beta_0\) and the moderator \(x\) multiplied by the coefficient \(\beta_1\). \(x\) is dummy-coded, where indicates that the observed \(y_i\) is an effect size estimate of compared to control, while indicates that \(y_i\) is an effect size estimate of the compared to control.

Tocilizumab vs. Sarilumab

Here are the priors for this comparison. Of note, the “Optimistic for Sarilumab” and “Optimistic for Tocilizumab” models were the ones originally preregistered as main models. Yet they are presented as supplementary in the peer-reviewed article of this project.

\[ \begin{align*} \beta_0 & \sim \operatorname{Normal}(0, 0.75^2) \tag{Priors}\\ \tau & \sim \operatorname{Log-Normal}(-1.975, 0.67^2)\\ \beta_{1[Skeptical]} & \sim \operatorname{Normal}(0, 0.354^2)\\ \beta_{1[Sarilumab]} & \sim \operatorname{Normal}(-0.049, 0.193^2)\\ \beta_{1[Tocilizumab]} & \sim \operatorname{Normal}(0.049, 0.193^2)\\ \beta_{1[Vague]} & \sim \operatorname{Normal}(0, 4^2) \end{align*} \]

## Means

# Figure S3 in https://doi.org/10.1101/2021.06.18.21259133;
# version posted June 22, 2021
reciprocal_remap_cap = 1/1.05
mean_log_sari = log(reciprocal_remap_cap)

mean_skeptical = log(1)

## SDs

# Assuming a mean = log(1), what is the SD that yields a distribution
# with 0.025 probability density below log(0.5)?
sari_sd_skeptical = (mean_skeptical - log(0.5))/qnorm(1-0.025)

# Assuming a mean = log(0.952), what is the SD that yields a distribution
# with 0.4 probability density above log(1)?
sari_sd_optimistic_sari = (log(1) - mean_log_sari)/qnorm(1-0.4)

# Assuming a mean = log(1/0.952) = -log(0.952), what is the SD that yields a
# distribution with 0.4 probability density below log(1)?
sari_sd_optimistic_toci = (-mean_log_sari - log(1))/qnorm(1-0.4)

sari_sd_vague = 4

# General model characteristics 

## Formula
mf = 
  formula(yi | se(sei) ~ 0 + Intercept + treatment + (1 | study))

## Priors for beta_0 and tau

# tau:
# logmean
# -1.975

# logsd
# 0.67

priors = 
  brms::prior(normal(0, 0.75), class = "b", coef = "Intercept") +
  brms::prior(lognormal(-1.975, 0.67), class = "sd", group = "study")

Fit models

## Skeptical model

# sari_sd_skeptical
# 0.353653

m_sari_skeptical =
  brms::brm(
    
  data = d_logOR |> dplyr::filter(treatment != "baricitinib"),
  family = gaussian,
  formula = mf,
  
  prior = priors +
    brms::prior(normal(0, 0.353653), class = "b", coef = "treatmentsarilumab"),
  sample_prior = TRUE,
  
  control = list(adapt_delta = .96),
  backend = "cmdstanr", # faster
  warmup = 2000,
  iter = 4000,
  seed = 123,
  cores = parallel::detectCores(),
  chains = 4,
  
  file = here::here("output", "fits", "mr_sari_skeptical.Rds"),
  file_refit = "on_change"
)

## "Optimistic for Sarilumab" model

# mean_log_sari
# -0.04879016

# sari_sd_optimistic_sari
# 0.1925823

m_sari_optimistic_sari =
  brms::brm(
    
  data = d_logOR |> dplyr::filter(treatment != "baricitinib"),
  family = gaussian,
  formula = mf,
  
  prior = priors +
    brms::prior(normal(-0.04879016, 0.1925823), class = "b", coef = "treatmentsarilumab"),
  sample_prior = TRUE,
  
  control = list(adapt_delta = .96),
  backend = "cmdstanr", # faster
  warmup = 2000,
  iter = 4000,
  seed = 123,
  cores = parallel::detectCores(),
  chains = 4,
  
  file = here::here("output", "fits", "mr_sari_optimistic_sari.Rds"),
  file_refit = "on_change"
)

## "Optimistic for Tocilizumab" model

# -mean_log_sari
# 0.04879016

# sari_sd_optimistic_toci
# 0.1925823

m_sari_optimistic_toci =
  brms::brm(
    
  data = d_logOR |> dplyr::filter(treatment != "baricitinib"),
  family = gaussian,
  formula = mf,
  
  prior = priors +
    brms::prior(normal(0.04879016, 0.1925823), class = "b", coef = "treatmentsarilumab"),
  sample_prior = TRUE,
  
  control = list(adapt_delta = .96),
  backend = "cmdstanr", # faster
  warmup = 2000,
  iter = 4000,
  seed = 123,
  cores = parallel::detectCores(),
  chains = 4,
  
  file = here::here("output", "fits", "mr_sari_optimistic_toci.Rds"),
  file_refit = "on_change"
)

## Vague model

m_sari_vague =
  brms::brm(
    
  data = d_logOR |> dplyr::filter(treatment != "baricitinib"),
  family = gaussian,
  formula = mf,
  
  prior = priors +
    brms::prior(normal(0, 4), class = "b", coef = "treatmentsarilumab"),
  sample_prior = TRUE,
  
  control = list(adapt_delta = .96),
  backend = "cmdstanr", # faster
  warmup = 2000,
  iter = 4000,
  seed = 123,
  cores = parallel::detectCores(),
  chains = 4,
  
  file = here::here("output", "fits", "mr_sari_vague.Rds"),
  file_refit = "on_change"
)

Tocilizumab vs. Baricitinib

Here are the priors for this comparison. Of note, the “Optimistic for Baricitinib” and “Optimistic for Tocilizumab” models were the ones originally preregistered as main models. Yet they are presented as supplementary in the peer-reviewed article of this project.

\[ \begin{align*} \beta_0 & \sim \operatorname{Normal}(0, 0.75^2) \tag{Priors}\\ \tau & \sim \operatorname{Log-Normal}(-1.975, 0.67^2)\\ \beta_{1[Skeptical]} & \sim \operatorname{Normal}(0, 0.354^2)\\ \beta_{1[Baricitinib]} & \sim \operatorname{Normal}(-0.105, 0.416^2)\\ \beta_{1[Tocilizumab]} & \sim \operatorname{Normal}(0.105, 0.416^2)\\ \beta_{1[Vague]} & \sim \operatorname{Normal}(0, 4^2) \end{align*} \]

## Means
mean_log_bari = log(0.9)

mean_skeptical = log(1)

## SDs

# Assuming a mean = log(1), what is the SD that yields a distribution
# with 0.025 probability density below log(0.5)?
bari_sd_skeptical = (mean_skeptical - log(0.5))/qnorm(1-0.025)

# Assuming a mean = log(0.9), what is the SD that yields a distribution
# with 0.4 probability density above log(1)?
bari_sd_optimistic_bari = (log(1) - mean_log_bari)/qnorm(1-0.4)

# Assuming a mean = log(1/0.9) = -log(0.9), what is the SD that yields a
# distribution with 0.4 probability density below log(1)?
bari_sd_optimistic_toci = (-mean_log_bari - log(1))/qnorm(1-0.4)

bari_sd_vague = 4

## Skeptical model

# bari_sd_skeptical
# 0.353653

m_bari_skeptical =
  brms::brm(
    
  data = d_logOR |> dplyr::filter(treatment != "sarilumab", study != "RECOVERY Bari"),
  family = gaussian,
  formula = mf,
  
  prior = priors +
    brms::prior(normal(0, 0.353653), class = "b", coef = "treatmentbaricitinib"),
  sample_prior = TRUE,
  
  control = list(adapt_delta = .96),
  backend = "cmdstanr", # faster
  warmup = 2000,
  iter = 4000,
  seed = 123,
  cores = parallel::detectCores(),
  chains = 4,
  
  file = here::here("output", "fits", "mr_bari_skeptical.Rds"),
  file_refit = "on_change"
)


## "Optimistic for Baricitinib" model

# mean_log_bari
# -0.1053605

# bari_sd_optimistic_bari
# 0.4158742

m_bari_optimistic_bari =
  brms::brm(
    
  data = d_logOR |> dplyr::filter(treatment != "sarilumab", study != "RECOVERY Bari"),
  family = gaussian,
  formula = mf,
  
  prior = priors +
    brms::prior(normal(-0.1053605, 0.4158742), class = "b", coef = "treatmentbaricitinib"),
  sample_prior = TRUE,
  
  control = list(adapt_delta = .96),
  backend = "cmdstanr", # faster
  warmup = 2000,
  iter = 4000,
  seed = 123,
  cores = parallel::detectCores(),
  chains = 4,
  
  file = here::here("output", "fits", "mr_bari_optimistic_bari.Rds"),
  file_refit = "on_change"
)

## "Optimistic for Tocilizumab" model

# -mean_log_bari
# 0.1053605

# bari_sd_optimistic_toci
# 0.4158742

m_bari_optimistic_toci =
  brms::brm(
    
  data = d_logOR |> dplyr::filter(treatment != "sarilumab", study != "RECOVERY Bari"),
  family = gaussian,
  formula = mf,
  
  prior = priors +
    brms::prior(normal(0.1053605, 0.4158742), class = "b", coef = "treatmentbaricitinib"),
  sample_prior = TRUE,
  
  control = list(adapt_delta = .96),
  backend = "cmdstanr", # faster
  warmup = 2000,
  iter = 4000,
  seed = 123,
  cores = parallel::detectCores(),
  chains = 4,
  
  file = here::here("output", "fits", "mr_bari_optimistic_toci.Rds"),
  file_refit = "on_change"
)

## Vague model

m_bari_vague =
  brms::brm(
    
  data = d_logOR |> dplyr::filter(treatment != "sarilumab", study != "RECOVERY Bari"),
  family = gaussian,
  formula = mf,
  
  prior = priors +
    brms::prior(normal(0, 4), class = "b", coef = "treatmentbaricitinib"),
  sample_prior = TRUE,
  
  control = list(adapt_delta = .96),
  backend = "cmdstanr", # faster
  warmup = 2000,
  iter = 4000,
  seed = 123,
  cores = parallel::detectCores(),
  chains = 4,
  
  file = here::here("output", "fits", "mr_bari_vague.Rds"),
  file_refit = "on_change"
)

Sensitivity analysis with all patients in RECOVERY Bari

Here we will use the post-hoc priors based on Karampitsakos et al.’s results for the “Optimistic for Baricitinib” and “Optimistic for Tocilizumab” models.

# Data provided by Karampitsakos et al. through email 

# Calculate crude log odds ratio

bari_toci_data = 
  metafor::escalc(
  measure = "OR", # log odds ratio
  
  # Baricitinib
  ai = 40,
  n1i = 125,
  
  # Tocilizumab
  ci = 50,
  n2i = 126
) %>%
  as_tibble() %>% 
  # Calculate standard error
  dplyr::mutate(sei = sqrt(vi))

## Mean
mean_log_bari = bari_toci_data$yi[1]

## SDs
bari_sd_optimistic_bari = bari_toci_data$sei[1]
bari_sd_optimistic_toci = bari_sd_optimistic_bari

## Skeptical model

# bari_sd_skeptical
# 0.353653

m_bari_skeptical_sens =
  brms::brm(
    
  data = d_logOR |> dplyr::filter(treatment != "sarilumab", study != "RECOVERY Bari (No Toci)"),
  family = gaussian,
  formula = mf,
  
  prior = priors +
    brms::prior(normal(0, 0.353653), class = "b", coef = "treatmentbaricitinib"),
  sample_prior = TRUE,
  
  control = list(adapt_delta = .96),
  backend = "cmdstanr", # faster
  warmup = 2000,
  iter = 4000,
  seed = 123,
  cores = parallel::detectCores(),
  chains = 4,
  
  file = here::here("output", "fits", "mr_bari_skeptical_sens.Rds"),
  file_refit = "on_change"
)


## "Optimistic for Baricitinib" model

# mean_log_bari
# -0.3350615

# bari_sd_optimistic_bari
# 0.2644288

m_bari_optimistic_bari_direct_sens =
  brms::brm(
    
  data = d_logOR |> dplyr::filter(treatment != "sarilumab", study != "RECOVERY Bari (No Toci)"),
  family = gaussian,
  formula = mf,
  
  prior = priors +
    brms::prior(normal(-0.3350615, 0.2644288), class = "b", coef = "treatmentbaricitinib"),
  sample_prior = TRUE,
  
  control = list(adapt_delta = .96),
  backend = "cmdstanr", # faster
  warmup = 2000,
  iter = 4000,
  seed = 123,
  cores = parallel::detectCores(),
  chains = 4,
  
  file = here::here("output", "fits", "mr_bari_optimistic_bari_direct_sens.Rds"),
  file_refit = "on_change"
)

## "Optimistic for Tocilizumab" model

# -mean_log_bari
# 0.3350615

# bari_sd_optimistic_toci
# 0.2644288

m_bari_optimistic_toci_direct_sens =
  brms::brm(
    
  data = d_logOR |> dplyr::filter(treatment != "sarilumab", study != "RECOVERY Bari (No Toci)"),
  family = gaussian,
  formula = mf,
  
  prior = priors +
    brms::prior(normal(0.3350615, 0.2644288), class = "b", coef = "treatmentbaricitinib"),
  sample_prior = TRUE,
  
  control = list(adapt_delta = .96),
  backend = "cmdstanr", # faster
  warmup = 2000,
  iter = 4000,
  seed = 123,
  cores = parallel::detectCores(),
  chains = 4,
  
  file = here::here("output", "fits", "mr_bari_optimistic_toci_direct_sens.Rds"),
  file_refit = "on_change"
)

## Vague model

m_bari_vague_sens =
  brms::brm(
    
  data = d_logOR |> dplyr::filter(treatment != "sarilumab", study != "RECOVERY Bari (No Toci)"),
  family = gaussian,
  formula = mf,
  
  prior = priors +
    brms::prior(normal(0, 4), class = "b", coef = "treatmentbaricitinib"),
  sample_prior = TRUE,
  
  control = list(adapt_delta = .96),
  backend = "cmdstanr", # faster
  warmup = 2000,
  iter = 4000,
  seed = 123,
  cores = parallel::detectCores(),
  chains = 4,
  
  file = here::here("output", "fits", "mr_bari_vague_sens.Rds"),
  file_refit = "on_change"
)

Model diagnostics

# Load function
source(here::here("functions", "diag_plot.R"))

## * The library is already synchronized with the lockfile.

Meta-analyses

Tocilizumab model

diag_plot(model = ma_toci,
          pars_list = c("b_Intercept", "sd_study__Intercept"),
          ncol_trace = 4)

brms::pp_check(ma_toci,
               type = "dens_overlay",
               ndraws = 50)

brms::pp_check(ma_toci,
               type = "ecdf_overlay",
               ndraws = 50)

Sarilumab model

diag_plot(model = ma_sari,
          pars_list = c("b_Intercept", "sd_study__Intercept"),
          ncol_trace = 4)

brms::pp_check(ma_sari,
               type = "dens_overlay",
               ndraws = 50)

brms::pp_check(ma_sari,
               type = "ecdf_overlay",
               ndraws = 50)

Baricitinib (sensitivity) models

diag_plot(model = ma_bari,
          pars_list = c("b_Intercept", "sd_study__Intercept"),
          ncol_trace = 4)

brms::pp_check(ma_bari,
               type = "dens_overlay",
               ndraws = 50) + coord_cartesian(y = c(0, 5))

brms::pp_check(ma_bari,
               type = "ecdf_overlay",
               ndraws = 50)

Baricitinib models

diag_plot(model = ma_bari_sens,
          pars_list = c("b_Intercept", "sd_study__Intercept"),
          ncol_trace = 4)

brms::pp_check(ma_bari_sens,
               type = "dens_overlay",
               ndraws = 50) + coord_cartesian(y = c(0, 5))

brms::pp_check(ma_bari_sens,
               type = "ecdf_overlay",
               ndraws = 50)

Meta-regressions

Tocilizumab vs. Sarilumab

“Skeptical” model

diag_plot(model = m_sari_skeptical,
          pars_list = c("b_Intercept", "b_treatmentsarilumab", "sd_study__Intercept"),
          ncol_trace = 4)

brms::pp_check(m_sari_skeptical,
               type = "dens_overlay_grouped",
               group = "treatment",
               ndraws = 50)

brms::pp_check(m_sari_skeptical,
               type = "ecdf_overlay_grouped",
               group = "treatment",
               ndraws = 50)

“Optimistic for Sarilumab” model

diag_plot(model = m_sari_optimistic_sari,
          pars_list = c("b_Intercept", "b_treatmentsarilumab", "sd_study__Intercept"),
          ncol_trace = 4)

brms::pp_check(m_sari_optimistic_sari,
               type = "dens_overlay_grouped",
               group = "treatment",
               ndraws = 50)

brms::pp_check(m_sari_optimistic_sari,
               type = "ecdf_overlay_grouped",
               group = "treatment",
               ndraws = 50)

“Optimistic for Tocilizumab” model

diag_plot(model = m_sari_optimistic_toci,
          pars_list = c("b_Intercept", "b_treatmentsarilumab", "sd_study__Intercept"),
          ncol_trace = 4)

brms::pp_check(m_sari_optimistic_toci,
               type = "dens_overlay_grouped",
               group = "treatment",
               ndraws = 50)

brms::pp_check(m_sari_optimistic_toci,
               type = "ecdf_overlay_grouped",
               group = "treatment",
               ndraws = 50)

“Vague” model

diag_plot(model = m_sari_vague,
          pars_list = c("b_Intercept", "b_treatmentsarilumab", "sd_study__Intercept"),
          ncol_trace = 4)

brms::pp_check(m_sari_vague,
               type = "dens_overlay_grouped",
               group = "treatment",
               ndraws = 50)

brms::pp_check(m_sari_vague,
               type = "ecdf_overlay_grouped",
               group = "treatment",
               ndraws = 50)

Tocilizumab vs. Baricitinib

“Skeptical” models

diag_plot(model = m_bari_skeptical,
          pars_list = c("b_Intercept", "b_treatmentbaricitinib", "sd_study__Intercept"),
          ncol_trace = 4)

brms::pp_check(m_bari_skeptical,
               type = "dens_overlay_grouped",
               group = "treatment",
               ndraws = 50) + coord_cartesian(y = c(0, 4))

brms::pp_check(m_bari_skeptical,
               type = "ecdf_overlay_grouped",
               group = "treatment",
               ndraws = 50)

diag_plot(model = m_bari_skeptical_sens,
          pars_list = c("b_Intercept", "b_treatmentbaricitinib", "sd_study__Intercept"),
          ncol_trace = 4)

brms::pp_check(m_bari_skeptical_sens,
               type = "dens_overlay_grouped",
               group = "treatment",
               ndraws = 50) + coord_cartesian(y = c(0, 4))

brms::pp_check(m_bari_skeptical_sens,
               type = "ecdf_overlay_grouped",
               group = "treatment",
               ndraws = 50)

“Optimistic for Baricitinib” models

diag_plot(model = m_bari_optimistic_bari,
          pars_list = c("b_Intercept", "b_treatmentbaricitinib", "sd_study__Intercept"),
          ncol_trace = 4)

brms::pp_check(m_bari_optimistic_bari,
               type = "dens_overlay_grouped",
               group = "treatment",
               ndraws = 50) + coord_cartesian(y = c(0, 4))

brms::pp_check(m_bari_optimistic_bari,
               type = "ecdf_overlay_grouped",
               group = "treatment",
               ndraws = 50)

diag_plot(model = m_bari_optimistic_bari_direct_sens,
          pars_list = c("b_Intercept", "b_treatmentbaricitinib", "sd_study__Intercept"),
          ncol_trace = 4)

brms::pp_check(m_bari_optimistic_bari_direct_sens,
               type = "dens_overlay_grouped",
               group = "treatment",
               ndraws = 50) + coord_cartesian(y = c(0, 4))

brms::pp_check(m_bari_optimistic_bari_direct_sens,
               type = "ecdf_overlay_grouped",
               group = "treatment",
               ndraws = 50)

“Optimistic for Tocilizumab” models

diag_plot(model = m_bari_optimistic_toci,
          pars_list = c("b_Intercept", "b_treatmentbaricitinib", "sd_study__Intercept"),
          ncol_trace = 4)

brms::pp_check(m_bari_optimistic_toci,
               type = "dens_overlay_grouped",
               group = "treatment",
               ndraws = 50) + coord_cartesian(y = c(0, 4))

brms::pp_check(m_bari_optimistic_toci,
               type = "ecdf_overlay_grouped",
               group = "treatment",
               ndraws = 50)

“Optimistic for Tocilizumab” model

diag_plot(model = m_bari_optimistic_toci_direct_sens,
          pars_list = c("b_Intercept", "b_treatmentbaricitinib", "sd_study__Intercept"),
          ncol_trace = 4)

brms::pp_check(m_bari_optimistic_toci_direct_sens,
               type = "dens_overlay_grouped",
               group = "treatment",
               ndraws = 50) + coord_cartesian(y = c(0, 4))

brms::pp_check(m_bari_optimistic_toci_direct_sens,
               type = "ecdf_overlay_grouped",
               group = "treatment",
               ndraws = 50)

“Vague” models

diag_plot(model = m_bari_vague,
          pars_list = c("b_Intercept", "b_treatmentbaricitinib", "sd_study__Intercept"),
          ncol_trace = 4)

brms::pp_check(m_bari_vague,
               type = "dens_overlay_grouped",
               group = "treatment",
               ndraws = 50) + coord_cartesian(y = c(0, 4))

brms::pp_check(m_bari_vague,
               type = "ecdf_overlay_grouped",
               group = "treatment",
               ndraws = 50)

diag_plot(model = m_bari_vague_sens,
          pars_list = c("b_Intercept", "b_treatmentbaricitinib", "sd_study__Intercept"),
          ncol_trace = 4)

brms::pp_check(m_bari_vague_sens,
               type = "dens_overlay_grouped",
               group = "treatment",
               ndraws = 50) + coord_cartesian(y = c(0, 4))

brms::pp_check(m_bari_vague_sens,
               type = "ecdf_overlay_grouped",
               group = "treatment",
               ndraws = 50)